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Abstract — A direct simulation investigations of contaminant transport has been conducted. The simulations are based on 
the unsteady longitudinal advection diffusion contaminant decay equation , which is discretized using finite element 
Galerkin’s method with backward difference time formulation with the application of COMSOL Multiphysics Software 
Package. It was observed that the simulation results, which are illustrated pictorially via line graphs and contours, 
consistently represented the analytical solutions of Mebine and George (2011), thereby demonstrating interesting features of 
the problem effected by the variety of contaminant decay time -dependent sources. 
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I. Introduction 

Mebine and George (2011) investigated analytically the longitudinal advection-dispersion mathematical model with 
contaminant decay due to time-dependent sources. A variety of strategies for time -dependent sources were presented with the 
environmental objective to avoid intermittently high pollution levels downstream of the outfall. This was geared towards 
reducing environmental impact and minimising economic cost for downstream users of streams or groundwater polluted on 
purpose with non-conservative contaminants. 

The literature is rife with several studies of modalities, methods and advocates, emphasizing the importance of time- 
dependent release of contaminants and their dispersion, transport, migration processes or factors in order to gain insights on 
how to optimally control or reduce the impacts (Dhar and Sinha, 1989; Mazumder and Dalai, 2000; Baverman et al., 1995; 
Kolditz et al., 1998; Sandrin et al., 2001; Sen et al., 2002; Serrano, 2001; Mebine and Smith, 2006; Mebine and Smith, 2009). 
Rubin and Atkinson (2001) noted particularly that the spatial and temporal distribution of contamination sources have a 
significant impact on contamination transport and plume distribution in the subsurface. 

In recent times, the environmental, engineering and physical sciences span a broad range of issues in their disciplines such as 
open channel hydraulics, sediment transport, stratified flow phenomena, transport and mixing processes, water quality and 
atmospheric modelling, just to mention but a few. These topics are studied in various ways including theoretical analyses, 
physical model experiments, field studies, numerical modelling and numerical simulations. Today, the mathematical 
modelling and numerical solutions of some physical and engineering phenomena have been made easy due to the advent of 
less expensive but more powerful computers and sophisticated Commercial and Educational Software Packages with the 
development of user-friendly vis-a-vis easy-to-use rules and instructions yet robust modelling codes. Such packages include 
among others, MAPLE, MATHEMATIC A, and MATLAB, which are known for their capabilities for symbolic 
computations and numerical simulations; CORMIX Modelling Package ( http://www.cormix.info/) , which is a software 
system for the analysis, prediction, and design of marine outfall mixing zones resulting from a continuous point discharge of 
effluents into open coastal waters (Doneker and Jirka, 2007), and COMSOL Multiphysics Package, which has many vital 
Trade Mark statements such as “COMSOL Multiphysics™ provides a powerful platform for enhancing student knowledge 
by first confirming the analytical solutions based upon simplifying assumptions using numerical solutions, and then 
extending the problem to multi -dimensions with the addition of one or more Multiphysics; COMSOL Multiphysics™ 
provides a robust platform to solve complicated transport phenomena problems and to generate results that lead to a better 
understanding of the physical significance and their real-world applications; COMSOL Multiphysics™ application can 
generate a broad range of results needed for a better fundamental understanding of the problem being studied” (Vasilev et al. , 
2015). To this end, it is the objective of this paper to use COMSOL Multiphysics for the numerical simulations of 
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longitudinal advection dispersion contaminant decay mathematical model with time -dependent sources for the baseline 
comparisons of the analytical predictions in Mebine and George (2011). 

The sections that follow hereafter are respectively, brief discussion of some of the applications of COMSOL Multiphysics 
numerical simulations; mathematical formulations of the problem and dimensionless forms of the governing equations; 
source types and discussion of simulation results, and general concluding remarks of the results of the previous sections. 

II. Comsol Multiphysics Numerical Simulations 

COMSOL Multiphysics is a Commercial Software Package that uses a finite element solver with discretization by Galerkin’s 
Method for the purpose of stabilization to prevent spurious oscillations. The “extremely fine” settings have been used for the 
meshes with COMSOL backward difference formulation for the time setting. The COMSOL Multiphysics has been applied 
in the numerical simulations of many environmental and engineering problems including Transport phenomena (Vasilev et 
al., 2015); Navier-Stokes equations in supersonic flow over a double wedge airfoil (Kolluru and Gopal, 2012); Burning 
Candle (Gritter, 2010); Newtonian and Non-Newtonian Blood Flow (Siebert and Fodor, 2009); Topology Optimization of 
Heat Transfer and Fluid Flow Systems (Dede, 2009); Single Particle in Laminar Flow Regime of a Newtonian Liquid 
(Dinesh, 2009); Multiphase Fluid Flow Model in Porous Media (Diaz-Viera et al., 2008), and Velocity, Pressure and 
Temperature Distributions near a Stagnation Point in Planar Laminar Viscous Incompressible Flow (Kaufman and Gutierrez- 
Miravete, 2008). There is no gainsaying that these handful applications of COMSOL Multiphysics are pointers to its 
efficiencies, versatilities and capabilities in handling physically feasible problems. 

III. Mathematical Formulations 

The rate of contaminant transport in groundwater or river flow is governed by many factors including advection, dispersion, 
diffusion, dilution or mixing, retardation, and decay. These factors are usually incorporated into transport equations which 
describe transport phenomena such as heat transfer, mass transfer, fluid, waves, momentum transfer and many others. A 
generic transport equation that describes transport phenomena as in Siddig (2016) is stated as 

^ + V ■ f(t,x,</>,V</>)= g(t,x,</>), (1) 

dt 


where (/>,f,t,g and xare respectively the quantity to be transferred, flux, time, source and the Cartesian coordinates 

(x, y,z) and V = i — + j — + k—, and where i,j,k are unit vectors. Equation (1) is indicative that all transfer processes 

dx dy dz 

constitute certain conservation principles, flux expressions, constitutive equations, initial conditions, and boundary 
conditions to a particular environmental and engineering system geometry for all dependent variables in relation to the 
independent variables involved. This implies that the specific application of equation (1) depends on the issue of interest 
under investigation. 

The factors that characterize the transport of species are known for their specific roles such that the generic transport 
equation (1) could be simplified accordingly bearing in mind the necessary and sufficient transport physics. To this end, as in 
Mebine and George (2011), for consideration of one-dimensional domain contaminant transport or migration undergoing 
advection, dispersion, sorption isotherms and first-order contaminant decay, the governing equation is written as 


dc T7 dc d 2 c 

— + V — = D — - - Ac , 
dt dx dx 2 


( 2 ) 


where V = —,D = —,R = 1 + ^L 

r r e 


. The appropriate initial and boundary conditions (Israel-Cookey et al., 2002; Williams 


and Tomasko, 2008) are written as 


t < 0 : c(x, t) - 0, 

<3) 

[c(co, t) — > 0. 
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Here c(x,t) represents contaminant concentration; jc , longitudinal or downstream coordinate; t , time; U , advective 
velocity; k , dispersion coefficient; A , first-order contaminant decay, degradation, or kinetic reaction parameter which 
accounts for the overall water-bed exchange of the contaminant; R , accounts for contaminant retardation (equilibrium 
absorption) due to sorption and other physical processes; its values are generally less than 10; p b , bulk density; k d , 
contaminant distribution coefficient. The distribution coefficient is a physically meaningful parameter as it is a simple ratio 
of the amount sorbed to the amount dissolved in phase at equilibrium, and 6 , effective porosity. It is important to note that 
the cross-sectional averaged contaminant concentration c(v,?)is released at a rate c 0 g(^)from the point source at x = 0 for 
all t ; c 0 is the solute concentration after mixing over the cross-section of the stream; g(t ) is time-dependent function, which 
accounts for the time-dependent source. 

It is in order to render the equations (2, 3) dimensionless through the variables and parameters 



Vt c „ VL A 

— , C = — , Pe = — , A = 
L c n D 


LA 

V 


(4) 


Such that the equations (2, 3) are written as 

dC dC _J_ pC_ _ t 
dr dX Pe dX 2 

r <0:C(X,r) = 0, 

r>0:{ C J 0 ' r) -* (r) ' 

[c(oo,r)— >0, 


(5) 

( 6 ) 


where I is a characteristic length of interest; X,the reduced or dimensionless downstream distance; A, strength of 
contaminant decay, and Pe , the column Peclet number (van Genuchten and Alves, 1982). Physically, the Peclet number 
measures the characteristic time of the diffusion process (l 2 /d) relative to the convection process (L/v) . Therefore, while 
large Pe values represent less mixing at the ad vecti ve/convective front causing a sharper front, smaller Pe values represent 
more mixing causing a broader front. 

The equations (5, 6) are physically appropriate for investigating the effect of influent boundary, which varies over time due 
to variable discharge rate when considering the contaminant release and discharge problem in coastal industrial sites. This 
situation is realistic in the sense that generally industries or waste sites release pollution in a finite time period, either because 
industrial firms have a finite life or the pollution problem is controlled after a certain time with the awareness of the 
contamination or government regulation. Mebine and George (2011) utilized equations (5, 6) in proffering various analytical 
solutions, which could serve as toolkits for the testing of numerical experimentations and simulations. This is the definite 
objective of this paper with the application of COMSOL Multiphysics. 

To be at par with the results of Mebine and George (2011), the parametric values used for the simulations are 
t = 0.25, 0.50, 0.75, 1.00; Pe = 1000; A = 0.25, 0.50, 1.00, but for the longitudinal distance, X , it ranges as indicated on the line 
graphs. 

IV. CONT AMINAT SOURCE TYPES AND SIMULATION RESULTS 

The following time-dependent contaminant source types are considered for the numerical simulations: (i) accelerated source, 
(ii) continuous source, and (iii) exponential decay source. 

1.1 Accelerated Source 

The accelerated source has application in MHD free -convection and mass transfer flow on a porous medium in a rotating 
fluid by enhancing radiative heat transfer (Israel-Cookey et al., 2002). This is represented by g(r)=r. 
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Figure 1: Accelerated Source Line graphs for Concentration vs. Longitudinal Distance 
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Figure 2: Accelerated Source Surface Concentrations 
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1.2 Continuous (Impulsive) Source 

The continuous source depicts a known or constant injection concentration, where g(r) = l. In Laser heat source 
applications, Al-Khairy and Al-Ofey (2009) regarded the continuous source as a model of long duration of laser pulse when 
short times (of the order of few or tens of relaxation time) are considered. 
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Figure 3: Continuous Source Line graphs for Concentration vs. Longitudinal Distance 
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Figure 4: Continuous Source Surface Concentrations 
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1.3 Exponential Decay Source 

For exponentially decaying source, g(r)~ exp(-/?r) , where j3 is the source term decay constant which is different from the 
contaminant decay constant. Applications of exponentially decaying source are found in models of degradation of 
chlorinated compounds (Aziz et al., 2000). For radioactive waste disposal field and long-term behaviour of recalcitrant non- 
aqueous phase liquid (NAPL) spills, Williams and Tomasko (2008) reported that model physical processes are best 
conceptualized with exponentially decaying source. Another very important application of the exponentially decaying source 
is in the field of pharmaco -kinetics (Joshua, 2008) for the examination of drug absorption in the gastrointestinal track. 

Figures 1, 3 and 5 represents respectively the accelerated, continuous and exponential decay source line graphs of the 
simulated concentration distribution versus longitudinal distance at different pore volumes (time) r with Peclet number, 
Pe = 1000 and decay strength, A = 0.25, 0.50 and 1.00, respectively. The results excellently agree with those of the 
analytical results thereby replicating the characteristics of the accelerated, continuous and exponential decay sources, which 
were extensively discussed in Mebine and George (2011). However, it is observed that the line graphs of the continuous and 
exponential decay sources experienced some spurious or turbulent effects at the concentration peaks. This may be, 
physically, due to the cross dispersion coefficients that were used in the simulations. Equally, overshoots are visibly seen in 
the line graphs, which are reduced as the decay strength increases. From the Figures it is seen that the source term 
concentration is lower than the concentration behind the advective fronts. The introduction of the decay strength reduces the 
contaminant concentration behind the advective fronts with increasing longitudinal distance from the source. This is so 
because the contaminant at the advective front has undergone more degradation during transport. 

On the other hand, Figures 2, 4 and 6 indicates the accelerated, continuous and exponential decay sources surface 
concentration distributions for the decay strength, A = 0.25 for various pore volumes, t = 0.25, 0.50 and 1.00, respectively. 
It is observed that increasing pore volumes causes the contaminant concentrations to decrease everywhere in the flow (the red 
colour), thereby increasing the green colour (the safe zone). In effect, the downstream experienced minimum environmental 
impact of the non-conservative contaminant. 


V. Concluding Remarks 

The concentration distribution behaviour of non-conservative contaminants is depicted by the results of the simulations of the 
advection-diffusion-reaction mathematical model with time-dependent sources. The simulations have clearly shown that the 
choice of time-dependent contaminant release can have a persistent and marked effect upon the concentrations experienced 
far downstream, thereby making visible the characteristics and objectives of the problem. In addition, the simulations showed 
cross dispersions in terms of spurious or turbulent effects, which were not observed in the analytical solutions. Physically, the 
spurious or turbulent effects account for mixing. The overall testament of the investigation is that the simulations have, 
indeed, confirmed the fact that the analytical results of Mebine and George (2011) are usable for preliminary predictive tools 
for simulating contaminant migration in streams due to time -dependent contaminant releases, and they can attest to the 
effectiveness of COMSOL Multiphysics. 
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Figure 5: Exponential Source Line graphs for Concentration vs. Longitudinal Distance 
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